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Abstract. In this paper, we develop an estimator for the ver- 
tical flux of horizontal momentum with arbitrary beam point- 
ing, applicable to the case of arbitrary but fixed beam point- 
ing with systems such as the Poker Flat Incoherent Scatter 
Radar (PFISR). This method uses information from all avail- 
able beams to resolve the variances of the wind field in ad- 
dition to the vertical flux of both meridional and zonal mo- 
mentum, targeted for high-frequency wave motions. The es- 
timator utilises the full covariance of the distributed mea- 
surements, which provides a significant reduction in errors 
over the direct extension of previously developed techniques 
and allows for the calculation of an error covariance ma- 
trix of the estimated quantities. We find that for the PFISR 
experiment, we can construct an unbiased and robust esti- 
mator of the momentum flux if sufficient and proper beam 
orientations are chosen, which can in the future be opti- 
mized for the expected frequency distribution of momentum- 
containing scales. However, there is a potential trade-off be- 
tween biases and standard errors introduced with the new ap- 
proach, which must be taken into account when assessing the 
momentum fluxes. We apply the estimator to PFISR mea- 
surements on 23 April 2008 and 21 December 2007, from 
60-85 km altitude, and show expected results as compared 
to mean winds and in relation to the measured vertical veloc- 
ity variances. 

Keywords. Meteorology and atmospheric dynamics (Waves 
and tides; Instruments and techniques) - Radio science (Re- 
mote sensing) 


1 Introduction 

Geophysical vertical flux of horizontal momentum measure- 
ments were pioneered by Lhermitte (1983) and Vincent and 
Reid (1983), the former for studying the spectral variance 
in underwater tidal channels and the latter for assessing 
the vertical structure of wave motions in the atmosphere. 
In particular, Vincent and Reid (1983) developed a tech- 
nique for ground-based radars to estimate radial velocity 
variances and corresponding momentum fluxes in the meso- 
sphere and lower thermosphere (MLT) using dual co-planar 
narrow radar beams. This technique has since been widely 
used by HF and VHF radars (e.g., Reid and Vincent, 1987; 
Fritts and Vincent, 1987; Fukao et al., 1988; Reid et ah, 
1988; Fritts and Yuan, 1989; Fritts et al., 1990; Sato, 1990, 
1993, 1994; Tsuda et ah, 1990; Wang and Fritts, 1990, 1991; 
Hitchman et ah, 1992; Nakamura et ah, 1993; Murayama 
et ah, 1994; Murphy and Vincent, 1993, 1998). Similar ap- 
proaches have been extended to MF broad beam interfero- 
metric radars (Thorsen et ah, 1997), incoherent scatter radars 
(ISRs) operating at VHF and UHF (Fritts et ah, 1992; Fritts 
et ah, 2006; Janches et ah, 2006), and Doppler lidar systems 
(Acott, 2009). More recently. Hocking (2005) has general- 
ized the approach for application to meteor radars, a tech- 
nique that has been employed in several studies of gravity 
waves (GWs) (e.g., Antonita et ah, 2008; Clemesha et ah, 
2009; Placke et ah, 201 lb,a), validated by Fritts et ah (2010) 
in applications of the SAAMER meteor radar in studies of 
mean winds, tides and GWs, and used by Vincent et ah 
(2010) in an assessment of GW momentum flux measure- 
ments by meteor radars. In this paper, we outline a tech- 
nique for multi-beam radars to estimate the vertical flux of 
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Table 1 . Beam directions used for simulation and experiment. 



Beam 

map #1 

Beam 

map #2 

# 

Az. 

El. 

Az. 

El. 

1 

N/A 

90.0° 

N/A 

90.0° 

2 

4.23° 

73.73° 

0.0° 

85.0° 

3 

-1.81° 

79.57° 

45.0° 

85.0° 

4 

82.04° 

80.22° 

90.0° 

85.0° 

5 

91.92° 

74.59° 

135.0° 

85.0° 

6 

-88.82° 

80.73° 

180.0° 

85.0° 

7 

174.45° 

81.43° 

-90.0° 

85.0° 


horizontal momentum. While similar to the method of Hock- 
ing (2005), this technique is more applicable to the case of 
fixed-look directions determined a priori, rather than to the 
case of line-of-sight wind measurements from random me- 
teor trail locations. This technique is particularly applicable 
to phased array radars such as the Advanced Modular Inco- 
herent Scatter Radar (AMISR) class of radars. 


2 Momentum flux estimation - technique 

2.1 Measurements 

ISRs are able to measure the line-of-sight (LOS) or radial 
ion/neutral velocity within a small volume in the MLT. A 
component velocity measurement j, made by an arbitrary 
beam, can be written as: 

Vj = a,/ -v + ej (1) 

where v = (v e ,v n ,v z ) is the vector velocity in ge- 
ographic coordinates and aj = (cij, e , aj^ n , a.j tZ ) = 
(cos 0 j sin (p j , cos 9 j cos <p j , sh\9 j) is a simplified straight- 
line geometry vector with <pj and 9j being the azimuth 
(east of north) and elevation angles, respectively. (Note that 
throughout this paper, we will use the notation aj n , aj z 
to refer to the eastward, northward, and vertical components 
for a given beam j.) The measurement has an associated 
eiror, ej, with expected value (ej) — 0. 

Single instrument ISR experiments typically consist of 
measurements Vj distributed in altitude, space and/or time 
depending on the system, its steering dexterity, etc. We will 
consider cases applicable to the pulse-to-pulse beam steering 
of the Poker Flat Incoherent Scatter Radar (PFISR). 

2.2 Vector velocities 

As described by Heinselman and Nicolls (2008) and Nicolls 
et al. (2010), an arbitrary number of component velocity 
measurements can be used to estimate the velocity vector. 


This approach can be expressed in matrix form as 


Vl ’ 


d\,e ®\,n &2,z 


Ve 


e\ 

V 2 

— 

Cl 2,e Cl2,n d2 ,z 


Vn 

+ 

e 2 

V j. 


_ Cl j 7 e Cl j,n Cl _ 


_Vz_ 


J) _ 


( 2 ) 


ui os = Av + e. 


(3) 


A weighted linear least squares estimator for v is then, 
v = argmin(||ri os — ArH^J (4) 

V ' ' 


where 1 1 ■ | |w represents a weighted f 2 -norm with weight ma- 
trix given by W = C -1 , where C is the covariance matrix of 
the measurements. In this case, with uncorrelated errors, C 
is simply a diagonal matrix with diagonal elements given by 
(ej), the estimated variances of the measurements, where the 
notation (•} denotes an expected value, i.e., 


C = 


'(ej) 0 0 - 
_ 0 0 


(5) 


Solution to Eq. (4) is given by the normal equations (e.g., 
Aster et al., 2005; Tarantola, 2005), 


(A r C“ 1 A)v = A 7 'C“ 1 u l0S 


(6) 


or, rearranging, 

v — (A r C -1 A) -1 A r C -1 t>i os . (7) 


This solution is equivalent to the maximum a posteriori 
(MAP) estimator. If A 7 C -1 A is ill-conditioned, or if the ap- 
plication of a physical constraint is appropriate, a regular- 
ization/constraint term can be added. For example, this ap- 
proach has been applied to this problem by Heinselman and 
Nicolls (2008) and Butler et al. (2010). The error covariance 
matrix of the solution without a constraint term is given by 

Cv= (a^C^a)” 1 . (8) 

2.3 Velocity variances and momentum fluxes 


An estimator for the radial velocity perturbation is 

v'j = v j-( v j) (9) 

where we assume that after removal of the mean value, 
(vj) = 0. The variance of this quantity may be written as 

(v'f) = (vj a ) + (e)). (10) 
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This is a biased estimator (biased by (e 2 -)) of the variance. 
Similarly, the covariance of any two measurements (i ^ j ) 
can be written as, 


Wj) 


(Vi'vj 


( 11 ) 


where measurement errors have been assumed to be uncorre- 
lated. 

The quantity (i>- i/,} may be related to component fluxes: 

1 J 

(v'iv'j) = (v , l f)a itn aj Jl + (vf)ai, e aj t e + ( vf)a ir M jiZ ( 12 ) 
~i(^n^e) T Qj,e®i,n) T (^n^z) (dinCl j z T Uj in ^i,z) 

T(u e v^) (ai, e a j, z -\- a j,e®i,z) • 

For the case of i — j (variance measurement), 

(v/ 2 ) = (v' 2 )a 2 jn + (vf)a 2 je + (vf)a\ z + (13) 

2(v l nV l e )ajeaj'n+2(v l n v l z )ajj,aj' Z + 2(v l e v' z )aj, e aj' Z . 

The vertical fluxes of horizontal momentum correspond to 
the two last terms, {v' n v' z ) and (v' e v(), which we would like to 
estimate. Several techniques exist in estimating these quanti- 
ties, including choosing look directions such that unwanted 
terms in Eq. (13) cancel (such as the co-planar beam method 
of Vincent and Reid (1983)) or by formulating a matrix of 
scaled variances, applicable to arbitrary pointing/meteor trail 
drift measurements (Hocking, 2005). Note that, throughout 
this section, it is implicitly assumed that the time-averaged 
quantities are constant over the probing volume - an assump- 
tion of horizontal homogeneity. This is assumed to be true of 
the background wind fields (driven by tidal and other low- 
frequency motions) as well as the amplitudes of the wave 
fields. However, covariance measurements involve the esti- 
mation of the expected value of non-homogenous quantities 
(i.e., v', the velocity perturbation of a single measurement); 
thus, one must be cautious when including these terms using 
distributed measurements. This will be discussed in more de- 
tail in the simulation of realistic wave fields (Sect. 3.3). 

More generally, if the full covariance matrix of measure- 
ments is able to be formed, then this can be related to the 
covariance matrix of the wind field by a standard covariance 
matrix transformation. 


£ = ADA' 


(14) 


Here, for a number of measurements N meas (corresponding, 
for example, to the number of look directions or beams), T 


is the /V meas x /V meas covariance matrix of the measurements. 


(0f> (v[v' 2 ) ... (v[v'j) 
(v\v' 2 ) (vf) ... (v' 2 i Yj) 


(15) 


the observation matrix A has been defined by Eq. (2), and the 
wind covariance matrix is given by 


D = 


(v 2 ) (v e v n ) (v e v z ) 
{v e V n ) ( V 2 ) (v n V z ) 
( v e v z ) (v n v z ) ( V 2 ) 


(16) 


A solution to Eq. (14) can be found directly, 
d = (a^"^ Sana'a)” 1 . (17) 


While useful, this approach is more readily solvable by trans- 
forming the two-dimensional £ and D matrices into one- 
dimensional matrices of observational and estimated quan- 
tities. For simplicity and convenience, in our case we simply 
solve a linear least squares problem for 6 unknowns, (v' 2 ), 
( v ' e 2 ), ( v ' z 2 ), (v' n v' e ), {v' n v' z ), {v' e v’ z ). In this case, we may write. 


<fi?> 

v[v' 2 ) 


V v = 


Wit y 

(v?) 


{' v' 2 v'j ) 


= B 


K) 

(V 2 ) 

(l) e V/i ) 

(v e v z ) 
. i v n v z ) 


— BD„ 


(18) 


(if) 


defining D v to be the vector of variances and co-variances. 
For A meas measurements/beams/look directions, the observa- 
tion matrix T v will be size A meas (A meas + l)/2. The geom- 
etry matrix B is size AW^ ( A meas + l)/2 x 6, with elements 
given by the factors in Eq. (12), 
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(19) 


The inclusion of the covariance terms is a unique aspect of 
this approach and distinguishes it from previous techniques. 
If those terms were unable to be estimated or were ignored, 
than the approach would reduce to that of Hocking (2005), 
except that Hocking (2005) formed vectors of weighted vari- 
ances and a corresponding geometry matrix, such that the 
system is exactly determined. That approach seems appro- 
priate for meteor scattering for arbitrary locations and times 
over a large field-of-view. In the way we have formulated it, 
the problem may be over- or underdetermined and is an ap- 
propriate formulation for fixed experimental beam pointing, 
wherein the B matrix will remain fixed for any given exper- 
iment and where the covariance matrix of the measurements 
can be computed. We will show via simulation that including 
the covariance terms reduces the standard errors significantly 
in the presence of non-zero measurement uncertainties. 

The estimator using the same approach as in the previous 
section is 


D=(B 7 ’E- 1 B)- | B r E- 1 E„. 


( 20 ) 


A great advantage of this approach is that the a posteriori er- 
ror covariance matrix (errors on the wind variances and mo- 
mentum fluxes) can be readily estimated. This matrix is given 
by 


E^CBE^Xr 1 - (2D 

Here, T, vv is the covariance matrix of the measurements E,, , 
which can be readily computed. The general covariance of 
any entry in the covariance matrix is given by. 


Co wiv-v'j, v' m v ' n ) = {v iv'jv'J',) - {v'iv'j){v' m v' n ) 

= «,><w+< w<w 

+ (ViV' m )( e j,n) S j,n + (VjV' n )(ej m ) S j,‘ 

+ (v'jV'm)( e tn) S i,n + (Vj v' n ) (e 2 JSi,n 


(22) 


+ (eL,)(el, 


\X: 


, i.2 


where we have used the fourth moment theorem for multi- 
variate normally distributed variables to expand the higher 


order moments, and 5,-j is the Kronecker delta (1 for i — j, 
0 otherwise). For example, we see that if i = j = m = n, then 

Cov(i>; 2 , S' 2 ) = 2 (<u/ 2 > + (e 2 )) 2 ; (23) 

For O' = j) ^ (m =n), 

Cov(v- 2 , v' 2 ) — 2{vi'v/) 2 ; (24) 

For (i — m) ^ (j — n ), 

Cov(C-t)'-, v-v'j) = ({v/ 2 ) + {, ef }) ((v/ 2 > + {ejfj+iv/v/) 2 ; 

(25) 

For i ^ j yL m ^ n, 

Cov(v-Vj , v' m v' n ) = (viv m )(vjv n ) + {viv n )(vj v m ); (26) 

and for i — m, j 

Cov(t)-t)'., v-v'j = ([v/ 2 ) + (e 2 >) (vjv n ) + (viV„)(viVj). 

(27) 

From these expressions, the full covariance matrix of the 
measurements can be computed, which can be propagated to 
compute the variance of the estimator. Misestimation of the 
covariance matrix and violation of the linear least squares 
assumptions of normality can lead to biased estimates of the 
momentum fluxes. As a result, in practice, it may be better to 
use a diagonal prior covariance matrix and only use the es- 
timated covariance matrix for estimates of the errors on the 
results. 

As a final note, the bias in Eq. (10) can cause biases in 
the estimated quantities, particularly if the covariance terms 
are included as described here. While the biases on the vari- 
ance estimates can be removed if they are known, we have 
found that it may be better to estimate them as well. The ap- 
proach described above can be augmented to account for the 
biases. In this approach, the forward model is more correctly 
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Zonal Distance (km) Zonal Distance (km) 

Fig. 1 . PFISR beam position maps used in this paper, each with 7 look directions and evaluated at 80 km altitude. Left: A standard experi- 
mental beam map that has been used in several PFISR experiments. Right: A theoretical beam map with smaller off-zenith angles. Azimuth 
and elevation angles are listed in Table 1. 


described to include the bias terms and Eq. ( 14) is modified 
such that 

E = ADA r + C (28) 

where C is the diagonal measurement error covariance ma- 
trix given by Eq. (5). Equation (18) can then also be modified 
to add Al m eas additional unknowns (( e \ }, (e?)...) with a cor- 
responding modification in the B matrix. 

3 Momentum flux estimation - simulation 
3.1 Simulation description 

In this section, we use the principles developed in Sect. 2 
to perform a simulation of momentum flux extraction. For 
these simulations, we use two experiment geometries, as in- 
dicated in Table 1 and plotted at 80 km altitude in Fig. 1. 
Beam map #1 is the experiment geometry used by the D- 
region PFISR measurements of Nicolls et al. (2010). This 
experiment consists of seven look directions, which are listed 
and numbered in Table 1, and include a vertical look direc- 
tion. Beam map #2 is a theoretical beam map with the same 
number of positions (7), but with the maximum off-zenith 
angle reduced to 5°, for reasons that will become apparent in 
the simulation of realistic wave fields. 

Given these beam directions, we can directly compute the 
a posteriori covariance matrix of the estimator using Eq. (21), 
assuming a diagonal a priori covariance matrix (in general, 
this is not a good assumption, as the a priori covariance ma- 
trix can be highly non-diagonal given both the properties of 


the wind field and the fact that measurements are reused to 
compute the full covariance matrix; this exercise is merely 
to show the influence of the observation matrix). In Fig. 2, 
we utilise this calculation to investigate the importance of 
each beam from Table 1 (beam map #1) in the estimation 
procedure. In the future, a similar procedure could be used 
to optimize experiment pointing geometry. The output er- 
ror matrix of the estimated wind covariances is estimated 
after removing each beam sequentially from the estimation 
procedure. The fractional increase in the variance of each 
estimated quantity is indicated in Fig. 2 (black dots). The 
red dots show the same calculation, but for the procedure 
neglecting the measurement covariance terms that were in- 
cluded in the estimator in Sect. 2. 

In general, we see that ( v 2 e ) and {v\) are most sensitive 
to the removal of beams (5,6) and (2,7), respectively. This 
is not surprising, as these beams have the lowest elevation 
angles and are the most non-redundant look directions. The 
quantity ( vl ) is very well determined and not that sensitive to 
any beams. This is somewhat surprising, given the presence 
of a vertical beam, but the fact that all beams have fairly high 
elevation angles means that they are all sensitive to vertical 
motions (and, note that this calculation does not take into ac- 
count anything about the magnitude of the component vari- 
ances, expected to be much smaller in the vertical direction). 
The horizontal wind covariance term, (v e v n ), is the most 
poorly determined parameter and roughly equally sensitive to 
all beams other than the vertical beam. Note that the chosen 
experiment geometry pointed beams in roughly cardinal di- 
rections, which could be improved/optimized in future exper- 
iments with consideration about sensitivity to the dominant 


www.ann-geophys.net/30/945/2012/ 


Ann. Geophys., 30, 945-962, 2012 


950 


M. J. Nicolls et al.: PFISR momentum flux 


{vl) 



Beam Number 


2 

10 


0 

10 


-2 

10 


-4 

10 

1 2 3 4 5 6 7 

(v e Vz) 


2 

10 


0 

10 


-2 

10 


-4 

10 

Beam Number 



<<> 



(vl) 


10 


10 


10 


10 


10 


10 


10 


10 



1 2 3 4 5 6 7 

B earn N umber 


Fig. 2. Using beam map #1, estimate of the fractional increase in the variance of each estimated wind variance/covariance term after remov- 
ing individual beams (indicated by the x-axis) from the estimation procedure, for (black) procedure including covariance terms and (red) 
procedure neglecting covariance terms. This process assumes a diagonal measurement error covariance matrix. The absolute value of the 
variance (for a measurement error covariance matrix given by the identity matrix) using all beams is indicated in the upper left of each panel 
for the procedure including the covariance terms. The ratio of the absolute values using all beams is indicated by the horizontal red dashed 
line. This turns out to be an overestimate, due mostly to the covariance of the measurements. 


scales and frequencies. The momentum flux terms, (v e v z ) 
and (v n v z ) seem fairly well-determined and largely sensitive 
to individual beams pointed in the westward and southward 
directions, respectively, again due to the non-redundancy of 
these beams. 

Comparing the results, including the covariance terms out- 
lined in Sect. 2 (black) and neglecting those terms (red), we 
see that in general including the covariance terms leads to 
less sensitivity in removing individual beams. The absolute 
magnitude of the variances is much smaller for the estimates 
including the covariances - a factor of ~10 in root-mean- 
square error for the term (v e v n ), and a factor of ~ 1.5-3 for all 
other terms. However, these estimates neglect the role of non- 
diagonal measurement error covariance matrix, which will be 
included next when examining the results of simulations. 

To investigate the statistical properties and any biases as- 
sociated with momentum flux extraction using the technique 
developed, we performed two types of simulations: the first 
(Sect. 3.2), to investigate the statistical aspects of momentum 
flux extraction, and the second (Sect. 3.3) to simulate realis- 
tic wave fields. 


3.2 Simulation 1: statistical properties of momentum 
flux extraction 

For the purposes of an initial simulation designed to inves- 
tigate the statistical properties of momentum flux extraction, 
we use the beam configuration described in the previous sec- 
tion with an assumed, non-essential background wind envi- 
ronment v. We then perform the following steps for N i iter- 
ations of the simulation: 

- Generate a 3 x 3 matrix, L, with elements drawn from a 
uniform distribution from —5 to 5 m s' 1 . 

- Use this matrix to generate a symmetric positive- 
semidefinite covariance matrix L 7 L which now rep- 
resents the wind variances (( v 2 ), (u„) and (vl), di- 
agonal elements) and covariances ((v e v n ), (v e v z ) and 
(v n v z ), off-diagonal elements). This approach will gen- 
erate a wind covariance matrix with variance terms with 
a median value of ~25 m 2 s -2 and ranging up to ~50- 
70m 2 s' 2 . Covariance terms will have a median value 
close to 0 m 2 s , ranging from ~ —50 to 50 m 2 s~ 2 . 

- For N 2 iterations, generate a time series of wind vec- 
tors with length /V 3 points (corresponding to the number 
of samples used to estimate the variances, covariances 
and wind vectors), drawn from the multi-variate normal 
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Fig. 3. Simulation results. Histograms of the biases (true value subtracted from the mean result of the simulation over N 2 iterations) for each 
wind field variance/covariance measurement. Blue lines correspond to a e = 2 m s - 1 and red lines correspond to a e = 5 m s~ 1 . Solid lines 
correspond to the standard momentum flux approach, neglecting the covariance terms. Dashed lines follow the approach outlined in this 
paper including the covariance terms and subtracting the expected value of the measurement error variances from the estimated variances. 
Legend values show the median and standard deviation of each dataset. 


distribution with covariance matrix L J L and mean v. 
Note that the resulting uncertainties on the wind vari- 
ances and covariances will be highly dependent on Nj, 
as shown and discussed later. 

- For each look direction, generate observed line-of-sight 
wind components using Eq. (9) with errors drawn from 
the zero-mean normal distribution with variance ( e 2 ). 

- Utilise Eq. (20) to estimate the momentum flux and 
wind component variances, as well as estimated mea- 
surement errors (e 2 ). 

- Estimate the means and standard deviations of the vari- 
ance and covariance terms over the N 2 iterations, which 
can be used to compute the standard error on the esti- 
mate as well as any biases on the results. 

We have run simulations as described above using N 1 = 
N 2 = 1000 and /V 3 = 60, with a e = sj ( e 2 ) equal to 0, 2 and 


5 m s~ 1 . Results are shown in Fig. 3, where each panel shows 
a histogram of the difference between the true value and the 
estimated value (averaged over the N 2 iterations). This rep- 
resents a histogram of the bias of the estimator. Here, blue 
lines correspond to the cr e = 2 ms -1 case and red lines cor- 
respond to the o e — 5ms ' 1 case. Solid lines correspond to 
the standard momentum flux approach, neglecting the covari- 
ance terms, whereas dashed lines correspond to the technique 
described in this paper. 

The variances, (v 2 ), {u 2 } and (i> 2 ), are reasonably well de- 
termined for the small error case, especially in the vertical 
direction. As already mentioned, this is due to the high el- 
evation angles of the beams and the presence of a vertical 
beam. However, the variances are clearly biased (in all three 
directions), which is most clearly seen on examination of the 
(vp histogram where the red curves are offset by 25 m 2 s -2 
(and the blue curves by 4 m 2 s -2 ). As discussed in Sect. 2, 
there is an expected bias in the line-of-sight variance by (e 2 ) 
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Fig. 4. Similar to Fig. 3, except curves correspond to histograms of standard deviations of wind variance and momentum flux estimates, 
computed over the Nj iterations. 


that propagates to the wind variances. This bias is removed in 
the approach described (dashed lines), because an estimate of 
the measurement error variance was removed from the mea- 
surement variances. The bias, however, does not affect the 
covariance terms. For these terms, we see unbiased results 
in all cases, relatively poorly determined for the (v e v n ) term 
and reasonably well determined for the vertical flux results. 

The approach using the covariance terms (dashed lines) 
significantly outperforms the traditional approach, as seen by 
the standard deviation of the biases - the distribution of the 
biases is much narrower. This is especially true for the (v e v n ) 
term, which is due to the additional information provided by 
the covariances. 

The errors in the approach are quantified in more detail 
in Fig. 4, where we plot histograms of the standard devi- 
ations of the estimates in a similar format to Fig. 3. Stan- 
dard deviations were computed using the N 2 iterations and 
represent the error in the estimate. Note that we earlier de- 
scribed an approach to estimate errors from the covariance 
matrix of the measurements, and we have confirmed that the 
standard deviations plotted in Fig. 4 agree with these esti- 


mates. Errors increase significantly as a e increases, but are 
much lower, for all parameters, for the technique using the 
covariance terms. For the momentum flux terms, the errors 
decrease by about 25 % for the a e = 2 m s -1 case and about 
40% for the a e = 5ms -1 . 

The conclusion that we draw from these simulations is that 
with (1) appropriate beam geometry and (2) sufficiently pre- 
cise measurements, we can make unbiased estimates with 
reasonable errors of the vertical flux of horizontal momen- 
tum. These errors can quickly become untenable if these cri- 
teria are relaxed. Errors are significantly reduced if the full 
covariance matrix of the measurements can be computed, as 
is the case for PFISR measurements. 

Given that the momentum flux is a measure of the cor- 
relation between horizontal and vertical wind variations, its 
determination is statistical in nature: a single measurement, 
even with no instrument errors, has at least 100% uncer- 
tainties, determined by the geometric mean of the horizontal 

and vertical wind variances (e.g., J (v^){v~)) and the frac- 
tional variance embodied in the momentum flux term (e.g.. 
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Fig. 5. The top row corresponds to a e = 0 m s 1 , the middle row to a e = 2 m s 1 , and the bottom row to a e = 5 m s 1 . Left: Fractional 
error in horizontal momentum flux determination as a function of horizontal to vertical wind correlation. As the correlation approaches 1, 
the fractional error approaches V2/y/~N^ = 0.18 (solid horizontal line). Right: Absolute error in horizontal momentum flux as a function of 

J i v x)( v z) (where results for (v e v z ) and ( v n v z ) have been merged). Blue and green lines with errorbars show mean and lcr spread for the 

black and red points, respectively, and gray line corresponds to J {v^){v\) /N-$. In all panels, black dots are from the simulation neglecting 
the covariance terms, where red dots are from the approach described in the text. 


Kudeki and Franke, 1998; Thorsen et al., 2000; Vincent et al., 
2010). For the case of this simulation, we have assumed that 
Nj = 60 statistically independent points exist with which to 
compute the line-of-sight variances and corresponding wind 
variances and co-variances. In Fig. 5, this issue is demon- 
strated. 

Figure 5 (left) shows the fractional error in estimating the 
horizontal momentum flux as a function of the correlation 


between the horizontal and vertical wind components, for the 
case of no measurement errors. In this case, the approach de- 
scribed here including the covariance terms agrees exactly 
with the traditional approach (or, that of Flocking, 2005, for 
arbitrary pointing). When the correlation is near 1, the frac- 
tional error approaches *J2/ = 0.18. As the correlation 
decreases, the fractional error increases dramatically, demon- 
strating the importance of this quantity (e.g., see Kudeki and 
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Fig. 6. Biases in the momentum fluxes and vertical velocity variance as a function of horizontal wavelength for beam map #1 (see Table 1) 
(top row) and beam map #2 (bottom row). Black curves show the solutions with only the variance terms included in the estimator, and red 
terms show results including the full covariance matrix. Simulations are with a e = 2ms _1 . 


Franke, 1998). The right-hand panel shows the absolute er- 
ror for this simulation, where it is quite clear for any given 

value of wind variances, (v^) (v^) / Nj (gray line) is a lower 
bound on the error (note that v x here is either component of 
the wind field, v e or v n ). The true error can be larger than 
this, depending on the magnitude of the correlation between 
the horizontal and vertical fluctuations. The importance of 
independent samples should be clear, as the only way to re- 
duce the uncertainties on the momentum flux estimates is 
to incoherently average in time or space. Adding measure- 
ment errors (lower two rows of Fig. 5) clearly makes things 
worse, as one would expect. However, including the covari- 
ance terms significantly reduces the fractional and absolute 
errors - ~3CM-0 % in the a e = 2ms -1 and ~60-70 % in the 
a e = 5 ms -1 . 

3.3 Simulation 2: realistic wave field 

In this section, we extend the simulations to include a real- 
istic wave field, mainly to investigate the effects of realistic 
horizontal variations in the perturbed radial velocity fields, 
similar to the approach taken by others (e.g., Fritts et al., 
2010, 2012; Vincent et al., 2010). We simulate a case which 
is similar to Case 1 of Fritts et al. (2010) and Fritts et al. 
(2012), in which we have fixed mean, diurnal tide and semi- 


diurnal tidal amplitudes for the background winds in addition 
to imposed gravity wave variations. 

In our modified Case 1, we have two gravity waves, with 
fixed amplitudes as given in Table 2, propagating orthogo- 
nally to one another with a horizontal wavelength and pe- 
riod drawn from a uniform random distribution from —500 
to 500 km, and 6 to 45 min, respectively. These waves will 
lead to fixed wind variances and momentum fluxes, given in 
Table 2. For the simulation, we follow the following steps 
for N i iterations of the simulation to generate a realistic time 
series of measurements: 

- Draw a wavelength and period for each of the two waves 
from a uniform random distribution, as described. 

- For N 2 iterations, generate observed line-of-sight wind 
components using the background winds, horizontally 
varying gravity wave fields and measurement errors 
drawn from the zero-mean normal distribution with 
variance cr 2 . 

- Estimate the momentum flux and wind component vari- 
ances for each of those iterations. 

- Estimate the means and standard deviations of the vari- 
ance and covariance terms over the Ni iterations, which 


Ann. Geophys., 30, 945-962, 2012 


www.ann-geophys.net/30/945/2012/ 


M. J. Nicolls et al.: PFISR momentum flux 


955 


6.5 
6 

5.5 
5 

4.5 
4 

3.5 


2.5 


( V e V 

iiiiiiiiijiiillliii)iiiiiiiiiiiiiiili|iiiiiiniii 

;!lilillll|||||||||il 

i I' 


( V n V s) 


||lli'll||||||il|lli||||| | ||j||H!ii|||||||||iii| 


2.5 


1.5 


!!!!!!;!;!!!!■!!!!!!!!!!!!!!!!!;!!!!!!!!!!!!!!!!! 1 

iiiliii!liiiiii"iiiij|j!jj;||;ii!!iiiijliiiiii!i 


0 . 5 . 


<«4> 




iiiiiiiiiiiiiiii, 


it 




.i" 


- 500 - 400 - 300 - 200-100 0 100 200 300 400 500 - 500 - 400 - 300 - 200-100 0 100 200 300 400 500 '- 500 - 400 - 300 - 200-100 0 100 200 300 400 500 



X x (km) \y (km) A x (km) 


Fig. 7. Same as Fig. 6 for standard errors. 


Table 2. Background winds, gravity wave parameters, and corresponding wind variances and momentum fluxes for the simulation in Sect. 3.3. 


Background winds 
(ms -1 ) 

GW1 

(ms -1 , km. min) 

GW2 

(ms -1 , km, min) 

Variances 
(in 2 s -2 ) 

(Cm- V M ) = (20, 10) 

IT 

II 

© 

O 

(j\ 

(k\ v\ w') = (0, 20, 2) 

«u 2 ), <u 2 ), <u 2 )) = (50, 200, 14.5) 

(U D , V D ) = (10, 10) 

X x = - 500 to 500 , X v = oo 

X x = oo , X y = - 500 to 500 

(v e v„) =0 

(.Usd, V sd ) = (50. 50) 

r = 6 to 45 

r = 6 to 45 

({v e v z ), (v n v z )) = (25,20) 


can be used to compute the standard error on the esti- 
mate as well as any biases on the results. 

The time series of the measurements corresponds to a two- 
hour interval with samples every 3 min, similar to the mea- 
surements that will be presented later. The simulation is per- 
formed at a single altitude of 80 km. 

It is crucially important to accurately estimate the co- 
variance terms, for which a direct covariance calculation 
is not suitable because of scale-dependent propagation de- 
lays across the field-of-view. To do so, we employ a 
method wherein we take the Fourier transform of the cross- 
covariance function, or the cross-power spectral density, be- 
tween measurement i and measurement j. We high-pass fil- 
ter the cross-spectrum in the frequency domain to eliminate 
low frequency motions with periods larger than an hour, 
which we are not able to address using the types of mea- 
surements targeted with this technique. We then integrate the 


cross-spectral density to obtain the magnitude of the covari- 
ance. The sign of the covariance is taken from the zero-lag 
of the cross-covariance function. Short-wavelength, spatially 
aliased waves may thus contribute to biased results, hence 
our desire to simulate a breadth of wave scales. 

The results for Case 1 are shown in Figs. 6 and 7. Fig- 
ure 6 shows the bias in the momentum fluxes and vertical 
velocity variance (average estimate minus true value) for the 
solutions including the variance terms in the estimator (black 
points) and the solutions including the covariance terms as 
outlined in this paper. Data are binned in terms of horizon- 
tal wavelength of the waves, since this was found to be a 
strong controlling factor on the nature of the bias. The top 
row of Fig. 6 utilises beam map #1 as indicated in Table 1. 
At large wavelengths, there is a clear bias in both methodolo- 
gies, approaching ~5-10 m 2 s -2 in the momentum fluxes. At 
large horizontal wavelengths, the vertical velocity variance 
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Fig. 8. Winds and momentum fluxes from 13:00-24:00 UT on 23 April 2008. Top: Zonal wind (black, v x ) and vertical flux of horizontal 
momentum (red, (v e v z ), top axis). Winds are two-hour averages plotted every hour, separated by 25 m s~ 1 . Momentum fluxes are computed 
over the two-hour period and are separated by 5 m 2 s -2 . Middle: Same for the meridional winds and momentum flux. Bottom: Vertical ve- 
locity (black, v z ) and vertical velocity variance (red, {v z ) 2 , top axis). Vertical winds are separated by lms^ 1 and vertical velocity variances 
by 5 nr s~ 2 . 


is unbiased as a result of accounting for the biased line- 
of-sight velocity variances. At small-horizontal wavelengths 
(<~100km) the method including the covariance terms be- 
comes significantly oppositely biased. This is clearly an ef- 
fect of the chosen geometry and scale-dependent propagation 
delays across the field-of-view. The bottom row of Fig. 6 
shows the same results for beam map #2, which has small 
off-zenith angles. Both the magnitude of the bias as well as 
the wavelength at which biases begins to deviate are signifi- 
cantly reduced. The importance of chosen beam geometry on 


the accuracy of the technique seems very clear. It should also 
be emphasized that the nature of any biases likely depends 
strongly on the wave field present at any given time. 

Figure 7 shows a similar plot for the standard errors on 
the estimates and demonstrates a tradeoff between accuracy 
and precision. While eliminating the covariance terms from 
the approach results in a less-biased solution, standard er- 
rors are significantly reduced by retaining those terms. Simi- 
larly, choosing a beam pattern with smaller off-zenith angles 
increases the accuracy of the estimates, but also results in 
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Fig. 9. Average momentum flux results from 13:00-24:00 UT on 23 April 2008. Top: Zonal (left) and meridional (right) wind (black lines) 
with error bars corresponding to the standard deviation of wind variation over the 1 1-h period. Red lines correspond to the averaged vertical 
flux of horizontal momentum in the zonal (left) and meridional (right) directions, with error bars corresponding to the standard error on the 
mean for the 1 1-h period. Blue lines correspond to the solution including only the variance terms. Bottom: Fraction of variance embodied in 
the momentum flux (black lines), averaged over the 1 1-h period, in the zonal (left) and meridional (right) directions. Red lines correspond to 
errors on the momentum flux estimates determined from the fitting algorithm averaged over the 1 1-h period (solid), average minimum error 
estimate using the geometric mean of the horizontal and vertical wind variances (dashed), and standard deviation of the data over the entire 
time period (dotted). 


larger standard errors. In the cases where biases are smaller 
than expected errors (such as in the case of larger measure- 
ment errors), it would make sense to trade increased preci- 
sion for accuracy. In any case, for a fixed experimental beam 
geometry, both solutions could be computed and compared 
to obtain the most quantitative picture possible. 

4 Results 

4.1 Observations on 23 April 2008 

We have applied the technique described in this paper to 
PFISR data from 23 April 2008, which has been described 
and analysed in detail by Nicolls et al. (2010). A large- 


scale, large-amplitude inertia-gravity wave was observed 
over Poker Flat during this period. For our purposes here, 
we have averaged incoherent scatter radar spectra in ~3 min 
integrations, corresponding to a sufficiently short average to 
capture the expected frequency range of velocity variability, 
and long enough to ensure that samples are approximately 
statistically independent. Spectra were fit as described in 
Nicolls et al. (2010) for spectral parameters including the 
line-of-sight drift velocity, taken to be equal to the compo- 
nent of the neutral wind along the radar line-of-sight. Line- 
of-sight velocities were averaged over 4 km in altitude (and 
analysed every 1 km). The covariance matrix of the mea- 
surements were computed over a two-hour period, corre- 
sponding to 38 independent variance estimates (which can 
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Fig. 10. Same as Fig. 8 for the period 18:00-24:00 UT on 21 December 2007. Note that momentum fluxes and vertical velocity variance 
spacing is 50 m 2 s~ 2 . 


be equated to IV 3 in the simulations in Sect. 3.2). As shown 
in the simulations and by Kudeki and Franke (1998), we ex- 
pect the minimum error on the zonal and meridional momen- 
tum flux estimates over the two-hour periods to correspond 

to y(i)g)(u 2 )/38 and y(i> 2 )(t> 2 }/38, respectively. 

In Fig. 8 , we show the two-hour averaged winds and mo- 
mentum fluxes in the zonal (top) and meridional (middle) di- 
rections, plotted every hour. The short vertical wavelength 
wave activity is clearly evident in the zonal and meridional 
winds. The momentum flux is well-behaved and is roughly 
anti-correlated with the wave perturbations, increasing in 
amplitude between 75 and 85 km. In the lower panel, we have 
plotted the averaged vertical velocities, along with the verti- 
cal velocity variances. The two-hour averaged vertical ve- 
locities fluctuate up to about 0.25 ms -1 (but more typically 


quite a bit less than that), and the variances for this dataset 
are quite small. 

In Fig. 9, we show averaged results for the 11-h period 
from 13:00-24:00 UT on 23 April. The top panels of Fig. 9 
show the averaged winds and momentum fluxes for this 11 -h 
period, determined by averaging each 2 -h estimate of the 
wind vector and momentum fluxes. Errorbars on the winds 
(black lines) correspond to standard deviations over the 11 -h 
period and are meant to indicate the variability of the back- 
ground winds embodied in Fig. 8 . The winds are predom- 
inantly westward and northward over this period, increas- 
ing with altitude and peaking at an altitude of about 80 km. 
The momentum fluxes are indicated by red lines, with error- 
bars corresponding to standard errors on the mean using the 
2-h estimates over the 11-h period. The average momentum 
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Fig. II. Same as Fig. 9 for the period 18:00-24:00 UT on 21 December 2007. 


fluxes are very small. The zonal flux is near zero except near 
65-67 km and 75-80 km when it increases in the westward 
direction to an averaged value of 0.5-1 m 2 s~ 2 . The merid- 
ional flux is southward over most of the altitude range, with 
an averaged value of about —0.75 to — 1 m 2 s~ 2 . Blue curves 
show the solutions including only the variance terms and 


show higher variability (as expected), but rough agreement 
with the solutions including the covariance terms. 

In the lower panels, we analyse the variability and ex- 
pected errors of the measurements. The black lines show the 
correlation of the horizontal and vertical winds, which reach 
values of ~40 % of the geometric mean of the wind vari- 
ances. Using the fact that there were at least 38 independent 
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estimates of the variance in each 2-h average, the clashed red 
lines show the expected (minimum) root-mean-square error 
on the momentum flux, again averaged over the 1 1-h period. 
The solid lines show the averaged errors determined from 
the mean-square error of the least-squares fitting process us- 
ing the covariance matrix of the measurements. The dotted 
lines show the standard deviation of the measurements over 
the 13:00-24:00 UT period. The curves indicate expected er- 
rors on the 2-hr momentum flux estimates of 0-2 m 2 s -2 . The 
curves have very similar altitude trends, with the estimated 
errors from the covariance matrix of the measurements some- 
what underestimating the variability of the data and the ex- 
pected minimum errors. 

4.2 Observations on 21 December 2007 

In Figs. 10 and 11, we show results in the same format 
for 21 December 2007 for the period from 18:00-24:00 UT 
when signal-to-noise ratios were sufficient to determine 
winds and momentum fluxes. In this case, data were averaged 
for 2 min before fitting for the radar spectra, and the analy- 
sis was done with 6-km altitude averaging (and computed 
every 3 km) to ensure robust estimates. In this case, verti- 
cal velocity excursions were significantly larger than in the 
previous case, with vertical velocity variances approaching 
~ 10 m 2 s -2 for individual two-hour windows, but more typi- 
cally quite a bit less than that. The momentum fluxes are cor- 
respondingly considerably larger (peaking at ~25 m 2 s -2 ), 
with fluctuations that mimic the dynamic background wind 
variations. However, in this case, it is important to realise 
that the horizontal and vertical wind variations are suffi- 
ciently large that it may not be possible to robustly deter- 
mine the momentum fluxes for a given two-hour window. As 
shown in the lower panels of Fig. 11, the expected averaged 
minimum root-mean-square error on the momentum flux is 
near 5m 2 s -2 . The averaged results in Fig. 11 show that 
the average momentum flux for this 6-h period is between 
— 10-0 m 2 s -2 in the zonal direction and ~ ± 1 0 m 2 s 2 in the 
meridional direction (varying with altitude). These average 
results should be significant given the statistical uncertain- 
ties and demonstrate larger momentum fluxes on this day as- 
sociated with the enhanced vertical velocities and variances. 
The predominant flux appears to be southwestward, oppo- 
site the direction of the background winds. Comparisons to 
the variance-only measurements show generally good agree- 
ment in the meridional direction, but indeterminate results in 
the zonal direction. 


5 Conclusions 

We have demonstrated a technique to estimate the covariance 
matrix of the wind field, including wind variances and mo- 
mentum fluxes, using arbitrary, distributed measurements of 
line-of-sight winds. This approach is applicable especially to 


systems with flexible, but fixed pointing, such as phased ar- 
ray radars like the Poker Flat ISR, and is useful especially for 
estimates of higher-frequency and shorter lived wave inputs 
rather than long-term averages. 

The approach is simply a generalization of previous tech- 
niques to include covariance terms. We demonstrated via 
simulation that the approach produces robust estimates of the 
wind variances and covariances that outperform the direct ex- 
tension of the co-planar beam technique in the presence of 
finite measurement uncertainties because of the inclusion of 
covariance terms in the forward model. Accurate estimation 
of the covariance terms is crucial to the approach, and scale- 
dependent propagation delays over the field-of-view of the 
system must be accounted for when estimating them. Mis- 
estimation of these terms and issues associated with experi- 
ment geometry, may lead to biased estimates of the momen- 
tum fluxes and wind variances. However, while complicating 
the approach, including the covariance terms provides a sig- 
nificant reduction in standard errors that may outweigh any 
potential biasing by misestimation. The covariance terms can 
easily be excluded to provide a less biased, but likely more 
error-prone, estimate of the momentum fluxes and wind vari- 
ances. 

Future work can optimize the beam positions, spatial av- 
eraging and time integration to take advantage of the ex- 
pected range of momentum-containing scales. In addition, 
other sources of potential biases, such as beam-pointing ac- 
curacy, should be investigated. 
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